www.gusucode.com > 基于粒子滤波的故障检测,使用似然函数作为检测函数 > 基于粒子滤波的故障检测,使用似然函数作为检测函数/code/FDI based on SIR likelihood/ceshi9.m

    close all;
clear all;
n = 1:600;%sample steps
stdw = sqrt(10);

npar = 500;%particle number
N = length(n);%N denotes the steps of samples
H=50:10:200;%threshold
faultalarm=zeros(1,16);%for every threshold,there is a missalarm rario reponding to it
A=50;F=300;%A denotes the number of simulation,F denotes the number of time points where there is fault 
j=0;
%N=500
% Generate the state process and observations
for k=1:16  %threshold from 50,60,70,......200
   switch k
   case 1
       h=50;
   case 2
       h=60;
   case 3
       h=70;
   case 4
       h=80;
   case 5
       h=90;
   case 6
       h=100;
   case 7
       h=110;
   case 8
       h=120;
   case 9
       h=130;
   case 10
       h=140;
   case 11
       h=150;
   case 12
       h=160;
   case 13
       h=170;
   case 14
       h=180;
   case 15
       h=190;
   case 16
       h=200;
end
   
    faultN=0;       %for every threshold,missN denotes the number of missing  alarm time points
for s=1:50        %for every threshold,50 times simulations is hold on
x0 = 0.1; %the prior position of x
c0 = 1;
b0=25;
xpath = zeros(1,N);
xmean=0;xvariance=0.1;
ymean=0;yvariance=1;
xnoise=gauss(xmean,xvariance,N); %the noise of station transite
ynoise=gauss(ymean,yvariance,N); 

b=b0;
xpath(1) = x0/2 + b*x0/(1+x0^2) +8*cos(0) +xnoise(1) ;
for i=2:N
    if i<301
        b=b0;
    else 
        b=b0*5;
    end
xpath(i) = xpath(i-1)/2 + b*xpath(i-1)/(1+xpath(i-1)^2) +8*cos(1.2*(i-1)) + xnoise(i);
end
zpath = 1/20*(xpath.^2) + ynoise;
% Particle filter with resampling
w = ones(npar,1)/npar;
xprev = randn(npar, 1);
SParMat = zeros(npar, N);
SXParMat = zeros(npar, N);
sxparpath = zeros(1,N);
likelihood=zeros(N,1);
D=zeros(N,1);
for i=1:N
    faultnumber=0;
xnext = drawpar(xprev, stdw, i);
xs = (xnext.^2)/20;
w = w.*(1/sqrt(2*pi)*exp(-((zpath(i)-xs).^2)/2));
l=1/sqrt(2*pi)*exp(-((zpath(i)-xs).^2)/2);
Li=sum(l)/npar;%i时刻的所有粒子似然函数值均值
likelihood(i)=Li;
Di=0;
if i<20
    for j=1:i
        Di=Di+(-(log(likelihood(j))));
    end
else 
    for j=i-20+1:i
        Di=Di+(-(log(likelihood(j))));
    end
end
D(i)=Di;
if i<301
    if D(i)>h
        faultnumber=faultnumber+1;
    else 
        faultnumber=faultnumber;
    end 
end
w = w/sum(w);
SParMat(:,i) = w;
SXParMat(:,i) = xnext;
sxparpath(i) = w'*xnext;
[xprev, w] = impResample(xnext, w);
end
faultN=faultN+faultnumber;%50次仿真中系统故障时残差值小于阈值的时间点总个数

end
 faultalarm(k)=faultN/(A*F);%计算每个阈值对应的漏报率，故障漏报率A是仿真总次数F是一次仿真中有故障时的时间点总个数
end


figure;

plot(H,faultalarm);